#include <riscv_vector.h>
#include<math.h>
#include<stdio.h>

#if !defined(DOUBLE)
#define VSETVL(n) vsetvl_e32m4(n)
#define FLOAT_V_T vfloat32m4_t
#define VLEV_FLOAT vle32_v_f32m4
#define VLSEV_FLOAT vlse32_v_f32m4
#define VSEV_FLOAT vse32_v_f32m4
#define VSSEV_FLOAT vsse32_v_f32m4
#define VFMACCVF_FLOAT vfmacc_vf_f32m4
#define VFMULVF_FLOAT vfmul_vf_f32m4
#else
#define VSETVL(n) vsetvl_e64m4(n)
#define FLOAT_V_T vfloat64m4_t
#define VLEV_FLOAT vle64_v_f64m4
#define VLSEV_FLOAT vlse64_v_f64m4
#define VSEV_FLOAT vse64_v_f64m4
#define VSSEV_FLOAT vsse64_v_f64m4
#define VFMACCVF_FLOAT vfmacc_vf_f64m4
#define VFMULVF_FLOAT vfmul_vf_f64m4
#endif

template <typename T>
inline void axpy_riscv1(const T* x, T* y, int len, const T alpha) {
  unsigned int jjj, lll;
  jjj = lll = 0;

  const unsigned int RISCV_STEP_SIZE = VSETVL(len);
  lll = len & ~(RISCV_STEP_SIZE-1);
  for (jjj = 0; jjj < lll; jjj += RISCV_STEP_SIZE) {
    VSEV_FLOAT(
      y + jjj,
      VFMACCVF_FLOAT(VLEV_FLOAT(y + jjj, RISCV_STEP_SIZE),alpha,VLEV_FLOAT(x + jjj, RISCV_STEP_SIZE),RISCV_STEP_SIZE), RISCV_STEP_SIZE);
  }

  for (; jjj < len; jjj++) {
    y[jjj] += alpha * x[jjj];
  }
}


template <typename T>
inline void axpy_riscv2(const T* x, T* y, int len, const T alpha) {
  unsigned int jjj, lll;
  jjj = lll = 0;
  unsigned int RISCV_STEP_SIZE = VSETVL(len);
  if (RISCV_STEP_SIZE <= len/2) {
			for (lll = 0, jjj=0; lll < len/(2*RISCV_STEP_SIZE); lll++, jjj+=2*RISCV_STEP_SIZE) {
        VSEV_FLOAT(y + jjj,VFMACCVF_FLOAT(VLEV_FLOAT(y + jjj, RISCV_STEP_SIZE),alpha,VLEV_FLOAT(x + jjj, RISCV_STEP_SIZE),RISCV_STEP_SIZE), RISCV_STEP_SIZE);
        VSEV_FLOAT(y + jjj + RISCV_STEP_SIZE,VFMACCVF_FLOAT(VLEV_FLOAT(y + jjj + RISCV_STEP_SIZE, RISCV_STEP_SIZE),alpha,VLEV_FLOAT(x + jjj + RISCV_STEP_SIZE, RISCV_STEP_SIZE),RISCV_STEP_SIZE), RISCV_STEP_SIZE);
			}
		}
  	for (; jjj < len; ) {
			RISCV_STEP_SIZE = VSETVL(len - jjj);
      VSEV_FLOAT(y + jjj,VFMACCVF_FLOAT(VLEV_FLOAT(y + jjj, RISCV_STEP_SIZE),alpha,VLEV_FLOAT(x + jjj, RISCV_STEP_SIZE),RISCV_STEP_SIZE), RISCV_STEP_SIZE);
			jjj += RISCV_STEP_SIZE;
		}
}

template <typename T>
inline void axpy(const T* x, T* y, int len, const T alpha) {
  unsigned int jjj;
  for (jjj = 0; jjj < len; jjj ++) {
    y[jjj] += alpha * x[jjj];
  }
}


template <typename T>
inline void axpy_noadd_riscv(const T* x, T* y, int len, const T alpha) {
  unsigned int jjj, lll;
  jjj = lll = 0;
  const unsigned int RISCV_STEP_SIZE = VSETVL(len);
  lll = len & ~(RISCV_STEP_SIZE-1);
  for (jjj = 0; jjj < lll; jjj += RISCV_STEP_SIZE) {
    VSEV_FLOAT(
      y + jjj,
      VFMULVF_FLOAT(VLEV_FLOAT(x + jjj, RISCV_STEP_SIZE), alpha, RISCV_STEP_SIZE),  RISCV_STEP_SIZE);

  }
  for (; jjj < len; jjj++) {
    y[jjj] = alpha * x[jjj];
  }
}